SUBROUTINE MPA_MYCRDS  &
 & (                    &
 & nx,    ny,    ns,    &
 & nxsh,  nysh,  nssh,  &
 & nxsh2, nysh2, nssh2, &

 & dims,  coords,       &

 & nxi, nxe,            &
 & nyi, nye,            &
 & nsi, nse,            &
 & nx0i, nx1e,          &
 & ny0i, ny1e,          &
 & ns0i, ns1e,          &
 & nx2e, ny2e,          &
 & nxsi, nxse,          &
 & nysi, nyse,          &
 & nssi, nsse,          &
 & nxs2i, nxs2e,        &
 & nys2i, nys2e,        &
 & nss2i, nss2e,        &
 & nx0si, nx1se,        &
 & ny0si, ny1se,        &
 & ns0si, ns1se,        &
 & ns0i_xproc, ns1e_xproc, &
 & ns0i_yproc, ns1e_yproc &
 & )

! use MPI

 use MPI_VARIABLES, only : &
 & decomp

 implicit none
 
 include 'mpif.h'

! Input variables
 integer(4), intent(in) :: nx,    ny,    ns
 integer(4), intent(in) :: nxsh,  nysh,  nssh
 integer(4), intent(in) :: nxsh2, nysh2, nssh2
 integer(4), intent(in) :: dims(1:3)
 integer(4), intent(in) :: coords(1:3)

! Output variables
 integer(4), intent(out) :: nxi, nxe
 integer(4), intent(out) :: nyi, nye
 integer(4), intent(out) :: nsi, nse
 integer(4), intent(out) :: nx0i, nx1e
 integer(4), intent(out) :: ny0i, ny1e
 integer(4), intent(out) :: ns0i, ns1e
 integer(4), intent(out) :: nx2e, ny2e
 integer(4), intent(out) :: nxsi, nxse
 integer(4), intent(out) :: nysi, nyse
 integer(4), intent(out) :: nssi, nsse
 integer(4), intent(out) :: nxs2i, nxs2e
 integer(4), intent(out) :: nys2i, nys2e
 integer(4), intent(out) :: nss2i, nss2e
 integer(4), intent(out) :: nx0si, nx1se
 integer(4), intent(out) :: ny0si, ny1se
 integer(4), intent(out) :: ns0si, ns1se
 integer(4), intent(out) :: ns0i_xproc, ns1e_xproc
 integer(4), intent(out) :: ns0i_yproc, ns1e_yproc

! Local variables
 integer(4) :: nstx, nsty, nsts
 integer(4) :: i, j, k
 integer(4) :: surplus

 if (decomp == 1) then
 
   nstx = int(float(nx)/float(dims(1)))
   nsty = int(float(ny)/float(dims(2)))
   nsts = int(float(ns)/float(dims(3)))

!  (nxi:nxe,nyi:nye,nsi:nse) - 
!  the basic domain of the process

   nxi = coords(1) * nstx + 1
   if (coords(1) + 1 < dims(1)) then
     nxe = (coords(1) + 1)*nstx
   else
     nxe = nx
   endif
 
   nyi = coords(2) * nsty + 1
   if (coords(2) + 1 < dims(2)) then
     nye = (coords(2) + 1) * nsty
   else
     nye = ny
   endif
 
   nsi = coords(3)*nsts + 1
   if (coords(3) + 1 < dims(3)) then
     nse = (coords(3) + 1)*nsts
   else
     nse = ns
   endif

 elseif (decomp == 2) then
 
!  Two-dimensional decomposition with smaller subdomains
!  near boundaries

   nstx = int(float(nx)/float(dims(1))) + 1
   surplus = nstx*dims(1) - nx
   nxi = 1
   nxe = nstx - int(float(surplus)/2.)
   if (coords(1) > 0) then
     do i = 1, coords(1)
       nxi = nxe + 1
       nxe = nxe + nstx
     enddo
   endif
   if (coords(1) == dims(1) - 1) nxe = nx
   
   nsty = int(float(ny)/float(dims(2))) + 1
   surplus = nsty*dims(2) - ny
   nyi = 1
   nye = nsty - int(float(surplus)/2.)
   if (coords(2) > 0) then
     do i = 1, coords(2)
       nyi = nye + 1
       nye = nye + nsty
     enddo
   endif
   if (coords(2) == dims(2) - 1) nye = ny

   nsts = int(float(ns)/float(dims(3)))
   nsi = coords(3)*nsts + 1
   if (coords(3) + 1 < dims(3)) then
     nse = (coords(3) + 1)*nsts
   else
     nse = ns
   endif

 endif

! print*, 'My x domain size is ', nxe - nxi
! STOP

! write (*,*) 'My dimension is', (nxe-nxi+1)*(nye-nyi+1)*ns, coords(1:3)

! (nxsi:nxse,nysi:nyse,nssi:nsse) - 
! the basic extended domain of the process (basic domain + exchange zone)

 nxsi = nxi - nxsh
 nxse = nxe + nxsh
 nysi = nyi - nysh
 nyse = nye + nysh
 nssi = nsi - nssh
 nsse = nse + nssh

! (nx0i:nx1e,ny0i:ny1e,ns0i:ns1e) - 
! the basic domain of the process, 
! extended at the boundaries of the model domain

 nx0i = nxi 
 nx1e = nxe 
 ny0i = nyi 
 ny1e = nye 
 ns0i = nsi 
 ns1e = nse 
 if (coords(1) == 0) nx0i = 0
 if (coords(1) == dims(1)-1) nx1e = nx + 1
 if (coords(2) == 0) ny0i = 0
 if (coords(2) == dims(2)-1) ny1e = ny + 1
 if (coords(3) == 0) ns0i = 0
 if (coords(3) == dims(3)-1) ns1e = ns + 1

 nx2e = nxe
 ny2e = nye
 if (coords(1) == dims(1)-1) nx2e = nx + 2
 if (coords(2) == dims(2)-1) ny2e = ny + 2

 nx0si = nx0i - nxsh
 nx1se = nx1e + nxsh
 ny0si = ny0i - nysh
 ny1se = ny1e + nysh
 ns0si = ns0i - nssh
 ns1se = ns1e + nssh

 nxs2i = nxi - nxsh2
 nxs2e = nxe + nxsh2
 nys2i = nyi - nysh2
 nys2e = nye + nysh2
 nss2i = nsi - nssh2
 nss2e = nse + nssh2
     
 nsts = int(float(ns+2)/float(dims(1))) ! Assumed dims(1) = dims(2)
 ns0i_xproc = coords(1)*nsts
 if (coords(1) /= dims(1) - 1) then
   ns1e_xproc = (coords(1) + 1)*nsts - 1
 else
   ns1e_xproc = ns + 1
 endif

 nsts = int(float(ns+2)/float(dims(2))) ! Assumed dims(1) = dims(2)
 ns0i_yproc = coords(2)*nsts
 if (coords(2) /= dims(2) - 1) then
   ns1e_yproc = (coords(2) + 1)*nsts - 1
 else
   ns1e_yproc = ns + 1
 endif


 return
 END SUBROUTINE MPA_MYCRDS
